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Abstract 

A recently presented method for the study of evolving self-gravitating 
relativistic spheres is applied to the description of the evolution of rela- 
tivistic polytropes. Two different definitions of relativistic polytrope, 
yielding the same Newtonian limit, are considered. Some examples 
are explicitly worked out, describing collapsing polytropes and bring- 
ing out the differences between both types of polytropes. 

Key words: Relativistic polytropes; post-quasi-static regime 



1 Escuela de Ffsica, Facultad de Ciencias, Universidad Central de Venezuela, Caracas, 
Venezuela. 

2 Postal address: Apartado Postal 80793, Caracas 1080A, Venezuela. e-mail: 
laherrer a@t elcel . net . ve 

3 Grupo de Ffsica Teorica, Departamento de Fi'sica, Facultad de Ciencias, Universidad 
de los Andes, Merida, Venezuela, e-mail: wbarreto@ula.ve 



1 



1 Introduction 



The use of polytropic equations of state in the study of the stellar structure 
has long and a venerable history [1, 2, 3] (and references therein). Its great 
sucess stems, mainly, from its simplicity and from the fact that it can be 
used to describe a large number of different situations. 

It is therefore not surprising that a great deal of work has been devoted 
to the study of polytropes in the context of general relativity [4, 5, 6, 7] (and 
references therein). Nevertheless, since the Lane-Emden equation, which is 
the cornerstone in the study of polytropic spheres, is based on the assumption 
of hydrostatic equlibrium, almost all works done so far (to our knowledge) 
on polytropic equations of state, concern spheres in hydrostatic equilibrium 
(collapsing "Newtonian" polytropes with n — 3, have been considered by 
Goldreich and Weber [8]). 

However, during their evolution, self-gravitating objects may pass through 
phases of intense dynamical activity, with time scales of the order of magni- 
tude of (or even smaller than) the hydrostatic time scale, and for which the 
static (or the quasi-static) approximation is clearly not reliable (e.g. the col- 
lapse of very massive stars [9] and the quick collapse phase preceding neutron 
star formation [10] (and references therein). In these cases it is mandatory 
to take into account terms which describe departure from equilibrium. Ac- 
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cordingly, it is our purpose in this work to study the evolution of polytropes 
out of hydrostatic equilibrium. 

For doing so, we shall make use of an approach for modeling the evolution 
of self-gravitating spheres, which does not require full numerical integration 
of time dependent Einstein equations [11]. The motivation for this, is based 
on the following argument: It is true that numerical methods [12] (and ref- 
erences therein) are enabling researchers to investigate systems which are 
extremely difficult to handle analytically. In the case of General Relativity, 
numerical models have proved valuable for investigations of strong field sce- 
narios and have been crucial to reveal unexpected phenomena [13]. Even 
specific difficulties associated with numerical solutions of partial differential 
equations in presence of shocks are being overpassed [14]. By these days what 
seems to be the main limitation for numerical relativity is the computational 
demands for 3D evolution, prohibitive in some cases [15]. Nevertheless, it 
is obviously simpler (in general) to solve ordinary than partial differential 
equations and furthermore purely numerical solutions usually hinder to catch 
general, qualitative, aspects of the process. Instead, the proposed method, 
starting from any interior (analytical or numerical) static spherically sym- 
metric ("seed") solution to Einstein equations, leads to a system of ordinary 
differential equations for quantities evaluated at the boundary surface of the 
fluid distribution, whose solution (numerical), allows for modeling the dy- 
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namics of self-gravitating spheres, whose static limit is the original "seed" 
solution. 

The approach is based on the introduction of a set of conveniently defined 
"effective" variables, which are effective pressure and energy density, and an 
heuristic ansatzs on the latter [11], whose rationale and justification become 
intelligible within the context of the post-quasistatic appproximation defined 
below. In the quasistatic approximation (see below), the effective variables 
coincide with the corresponding physical variables (pressure and density) 
and therefore the method may be regarded as an iterative method with each 
consecutive step corresponding to a stronger departure from equilibrium. In 
this work we shall restrain ourselves to the post-quasistatic level (see section 
4 for details). 

The fluid distribution under consideration will be assumed to be dissipa- 
tive. Indeed, dissipation due to the emission of massless particles (photons 
and/or neutrinos) is a characteristic process in the evolution of massive stars. 
In fact, it seems that the only plausible mechanism to carry away the bulk of 
the binding energy of the collapsing star, leading to a neutron star or black 
hole is neutrino emission [16]. Consequently, in this paper, the matter dis- 
tribution forming the self gravitating object will be described as a dissipative 
fluid, which in the equilibrium regime satisfies a polytropic equation of state. 
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On the other hand, in the treatment of radiative transfer within stellar 
objects, two different approximations are usually adopted: difussion and 
streaming out. 

In the diffusion approximation, it is assumed that the energy flux of 
radiation (as that of thermal conduction) is proportional to the gradient of 
temperature. This assumption is in general very sensible, since the mean free 
path of particles responsibles for the propagation of energy in stellar interiors 
is in general very small as compared with the typical length of the object. 
Thus, for a main sequence star as the sun, the mean free path of photons 
at the centre, is of the order of 2 cm. Also, the mean free path of trapped 
neutrinos in compact cores of densities about 10 12 g. cm. ~ 3 becomes smaller 
than the size of the stellar core [17, 18]. Furthermore, the observational 
data collected from supernovae 1987A indicates that the regime of radiation 
transport prevailing during the emission process, is closer to the diffusion 
approximation than to the streaming out limit [19]. 

However in many other circumstances, the mean free path of particles 
transporting energy may be large enough as to justify the free streaming ap- 
proximation. In this work, for simplicity, we shall consider only the streaming 
out limit. 

As we shall see, in the relativistic regime, two (at least) different defini- 
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tions of polytrope are possible, yielding the same Newtonian limit. We shall 
consider both of them, as possible "seed" equations of state and we shall 
contrast the patterns of evolution obtained from each case. 

The plan of the paper is as follows. In Section 2 we define the conventions 
and give the field equations and expressions for the kinematical and physical 
variables we shall use, in noncomoving coordinates. In Section 3 we present 
a brief review of the properties of Newtonian polytropes and discuss two 
possible generalizations to the relativistic regime. A resume of the proposed 
approach is presented in Section 4. In Section 5 the method is applied to 
the case when the "seed" equation of state is a relativistic polytrope and 
some examples are explicitly worked out. Finally a discussion of results is 
presented in Section 6. 

2 Relevant Equations and Conventions 
2.1 The field equations 

We consider spherically symmetric distributions of collapsing fluid, undergo- 
ing dissipation in the form of free streaming radiation, bounded by a spherical 
surface E. 

The line element is given in Schwarzschild-like coordinates by 
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ds 2 = e u dt 2 - e x dr 2 - r 2 (d6 2 + sin^dcf) 2 ) , (1) 

where v(t,r) and \(t,r) are functions of their arguments. We number the 
coordinates: x° = t; x 1 = r; x 2 = 9; x 3 = 0. We use geometric units and 
therefore we have c — G—1. 

The metric (1) has to satisfy the Einstein field equations 



G v » = -SttT;, (2) 

which in our case read [20]: 



-8< = 4 + e " A (^-7)' (3) 



-8JTT, 1 =-\ + e- A ( 4 + - I • 



-8ttT 2 = -8ttT 3 3 = - e — (2A + A(A - z>)) 



(4) 



4 

+ tl Uu" + z/ 2 - AV + 2^^) , (5) 



-8ttT 01 = --, (6) 
r 
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where dots and primes stand for partial differentiation with respect to t and 
r, respectively. 

In order to give physical significance to the T£ components we apply the 
Bondi approach [20]. Thus, following Bondi, let us introduce purely locally 
Minkowski coordinates (r,x,y,z) 

dr = e u l 2 dt ; dx = e x ^ 2 dr ; dy = rdO ; dz = rsinOdcj). 

Then, denoting the Minkowski components of the energy tensor by a bar, we 
have 

f ° = T ° ; T\ = T\ ■ Tl = T| ; f| = T 3 3 ; f 01 = e -("+ A )/ 2 T 01 . 

Next, we suppose that when viewed by an observer moving relative to these 
coordinates with proper velocity uj in the radial direction, the physical con- 
tent of space consists of a fluid of energy density p, radial pressure P and 
unpolarized radiation of energy density e traveling in the radial direction. 
Thus, when viewed by this (comoving with the fluid) observer the covariant 
tensor in Minkowski coordinates is 

/ p + e -e \ 

-e P + e 

P ' 

V OPj 

Then a Lorentz transformation readily shows that 
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1 — UT 



' I (8) 



Tl = T| = f 2 2 = f| = -P, (9) 



T 01 = e^)/ 2 f 01 = _ (P + P V^ +A)/2 _ e (*+A)/2 c (10 ) 

1 — CO 2 



with 



e^eil±^. (11) 
(1-w) 1 J 

Note that the coordinate velocity in the (t,r,9,(f>) system, dr/dt, is related 
to uj by 



w = ^: e (A-o/2. (12) 

dt 

Feeding back (7-10) into (3-6), we get the field equations in the form 



P + Puj 2 1 f 1 A / 1 A'\ 1 

1 — ur 8n \ r l \r z r \ 
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P + pu 2 
1 - co 2 



8n I r 2 




(14) 




(P + P) 



(u+X)/2 + e (u+X)/2 e 



8nr 



A 



(16) 



Observe that if v and A are fully specified, then (13-16) becomes a system 
of algebraic equations for the physical variables p, P, uj and e. 

At the outside of the fluid distribution, the spacetime is that of Vaidya, given 



where u is a coordinate related to the retarded time, such that u = constant 
is (asymptotically) a null cone open to the future and 1Z is a null coordinate 
(9kk = 0)- It should be remarked, however, that strictly speaking, the 
radiation can be considered in radial free streaming only at radial infinity. 

The two coordinate systems (t, r, 9, 0) and (u, 71, 9, <fi) are related at the 
boundary surface and outside it by 



by 




(17) 




(18) 
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Ti = r. 



(19) 



In order to match smoothly the two metrics above on the boundary sur- 
face r = rs(t), we must require the continuity of the first and the second 
fundamental form across that surface. Then it follows [11] 



e vs = l- 2M/R S , (20) 
e -A E = j _ 2M/i? E . (21) 

[Pk = 0, (22) 

where, from now on, subscript £ indicates that the quantity is evaluated 
at the boundary surface £. Next, it will be useful to calculate the radial 
component of the conservation law 



T u % = 0. (23) 

where 

Tfa, = {p + P) u^u v - Pg^ + dJu (24) 

with 



e -u/2 ue -X/2 



V(l"^ 2 ) 1/2 ' (1 -LU 2 ) l/ ' 



m ,0,0), (25) 
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Z" = (e'^ 2 , e~ A / 2 , 0, 0) , (26) 

where denotes the four velocity of the fluid, and P is a null outgoing 
vector. 

After tedious but simple calculations we get 

(-8*7?)' = (Tl - if) + W (r; - Jj) + ILI (a + f - f ) , (27) 
which in the static case becomes 

P' = -^(p + P), (28) 
which is the well known the Tolman-Oppenheimer-Volkoff equation. 

3 Newtonian and Relativistic Polytrope 

Although Newtonian polytropes are well known and examined in detail in 
most classical books on stellar structure [3] , we found it worthwhile to present 
here the very basic facts about its theory. 
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3.1 The Newtonian case 



As mentioned before, the theory of polytropes is based on the assumption of 
hydrostatic equlibrium, therefore the two starting equations are (remember 
that we are using geometric units) 

dP d(p 

Tr = (29) 



and 



with and p denoting the Newtonian gravitational potential and the mass 
(baryonic) density respectively. 

Combining the two equation above with the polytropic equation of state 

P = Kpl = Kpl +1/n , (31) 
one obtains the well known Lane-Emden equation (for 7^1) 

with 

r = CMo, (33) 



A °~ K(n + 1)' { ] 



(n-l)/n 
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K = Po/pOc, 



(35) 



where subscript c indicates that the quantity is evaluated at the centre, and 
the following boundary conditions apply: 



The boundary surface of the sphere is defined by £ = £„, such that ipo(£ n ) = 0. 

As it is well known, bounded configurations exist only for n < 5 and 
analytical solution may be found for n — 0, 1 and 5. 

It is also worth remembering that the polytropic equation of state may 
be used to model two different type of situations, namely: 

• When the polytropic constant K is fixed and can be calculated from 
natural constants. This is the case of a completely degenerate gas in 
the non-relativistic (7 = 5/3; n = 3/2) and relativistic limit (7 = 4/3; 
n = 3). 

• When K is a free parameter as for example in the case of isothermal 
ideal gas or in a completely convective star. 



#0 



(£ = 0) = 0; 



= 0) = 1. 



14 



3.2 The relativist ic case 



When considering the polytropic equation of state within the context of gen- 
eral relativity, two distinct expressions are often considered. In order to avoid 
confussion we shall differentiate them from the begining. Thus, the following 
two cases may be contemplated. 

3.3 Case I 

In this case the original polytropic equation of state is conserved 

P = Kpl +1/n , (36) 
then it follows from the first law of thermodynamics that 

^)-f = ^- < 37 > 

where T denotes temperature, a is entropy per unit of proper volume and M 
is the particle density, such that 

po = NrriQ. (38) 
Then for an adiabatic process it follows 



d(fy + Pd{jf) = 0, (39) 



which together with (36) leads to 



dpo 
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with 

7 = 1 + l/n. 
If 7 1, (40) may be easily integrated to give 

p = Cp + P/( 1 -l). 

In the non-relativistic limit we should have p — > p , and therefore C 
Thus, the polytropic equation of state amounts to 

p = po + P/( 7 - 1). 
It is worth noticing that the familiar "barotropic" equation of state 

p = P/( 7 -l), 

is a particular case of (42) with C — 0. 

In the very special case 7 = 1, one obtains 

„ i _ rf(p/po) 

whose solution is 

p = Plogpo + poC, 
or, if puting C = 1 from the non-relativistic limit 

p = Plogp + p . 
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From now on we shall only consider the 7 ^ 1 case. 
Next, let us introduce the following variables 

ot = Pjpc, (48) 

r = £/A (49) 
A 2 = Aixp c /[a(n + 1)], (50) 

C = Po/Poc (51) 
v(0 = m(r)A 3 /(47i Pc ), (52) 
where the mass function, as usually is defined by 

e~ x = 1 - 2m/ r. (53) 

Then the Tolman-Oppenheimer-Volkoff equation (28) becomes 

a£ (1 — a) + (n + l)ct0o 
and from the definition of mass function and equation (13) in the static case, 
we have 

m! = Aixr 2 p (55) 



or 

dv 



= CVot 1 -na + na^ ). (56) 
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In the Newtonian limit (a — > 0), (54) and (56) become 



and 

which are equivalent to the classical Lane-Emden equation 
3.4 Case II 



Sometimes it is assumed that the relativistic polytrope is defined by 

P = Kp 1+1/n , 
instead of (36). Then introducing 

i> n = p/pc 

related to ipo by 

i> n — V'o (1 - na + amp ). 
The TOV equation becomes 

at, 1 + 
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and from (55) 

* = «V. (64) 

In the Newtonian limit (a — > 0), the Lane-Emden equation is also recovered 
in this case, as it should be. 

Obviously both equations of state differ each other, specially in the highly 
relativistic regime. This can be verified by inspection of figures 1 and 2. 

4 The method 

Let us now give a brief resume of the method we shall use to describe the 
evolution of the relativistic polytrope. However before doing so some general 
considerations will be necessary. 

4.1 Equilibrium and quasi— equilibrium 

The simplest situation, when dealing with self-gravitating spheres, is that of 
equilibrium (static case). In our notation that means that oj = e = 0, all time 
derivatives vanishes, and we obtain the generalized Tolman-Oppenheimer- 
Volkoff equation (28). 

Next, we have the quasistatic regime. By this we mean that the sphere 
changes slowly, on a time scale that is very long compared to the typical time 
in which the sphere reacts to a slight perturbation of hydrostatic equilibrium, 
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this typical time scale is called hydrostatic time scale [3] (sometimes this time 
scale is also referred to as dynamical time scale, e.g. see [21]). Thus, in this 
regime the system is always very close to hydrostatic equilibrium and its 
evolution may be regarded as a sequence of static models linked by (16). 
This assumption is very sensible because the hydrostatic time scale is very 
small for many phases of the life of the star. It is of the order of 27 minutes 
for the Sun, 4.5 seconds for a white dwarf and 1CT 4 seconds for a neutron 
star of one solar mass and 10 Km radius. It is well known that any of the 
stellar configurations mentioned above, generally, change on a time scale that 
is very long compared to their respective hydrostatic time scales. 

However, as already mentioned, in some important cases, this approx- 
imation is not longer reliable, and one needs to consider departures from 
quasi-equilibrium. We shall describe such departures, in the post-quasi- 
static approximation defined below. 

4.2 The effective variables and the post— quasistatic ap- 
proximation 

Let us now define the following effective variables: 



= 1? = + (65) 
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p = - 




P + pu 2 
l-u 2 



+ €. 



(66) 



In the quasistatic regime the effective variables satisfy the same equation 
(28) as the corresponding physical variables (taking into account the contri- 
bution of e to the "total" energy density and radial pressure, whenever the 
free streaming approximation is being used). Therefore in the quasistatic 
situation (and obviously in the static too), effective and physical variables 
share the same radial dependence. Next, feeding back (65) and (66) into (13) 
and (14), these two equations may be formally integrated, to obtain: 



From where it is obvious that for a given radial dependence of the effec- 
tive variables, the radial dependence of metric functions becomes completely 
determined. 

With this last comment in mind, we shall define the post-quasistatic 
regime as that corresponding to a system out of equilibrium (or quasiequilib- 
rium) but whose effective variables share the same radial dependence as the 




(67) 




(68) 
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corrresponding physical variables in the state of equilibrium (or quasiequilib- 
rium). Alternatively it may be said that the system in the post-quasistatic 
regime is characterized by metric functions whose radial dependence is the 
same as the metric functions corresponding to the static (quasistatic) regime. 
The rationale behind this definition is not difficult to grasp: we look for a 
regime which although out of equilibrium, represents the closest possible sit- 
uation to a quasistatic evolution (see more on this point in the last Section). 

4.3 The algorithm 

Let us now outline the approach that we shall use: 

1. Take an interior solution to Einstein equations, representing a fluid 
distribution of matter in equilibrium, with a given 

Pst = p(r); P st = P(r). 

This static solution will be obtained in this work by integration of the 
relativistic Lane-Emden equations (54), (56) or (63), (64). 

2. Assume that the r dependence of P and p is the same as that of P st 
and p st , respectively. 

3. Using equations (67) and (68), with the r dependence of P and p, one 
gets m and v up to some functions of t, which will be specified below. 
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4. For these functions of t one has three ordinary differential equations 
(hereafter referred to as surface equations), namely: 

(a) Equation (12) evaluated on r = r s . 

(b) The equation relating the total mass loss rate with the energy flux 
through the boundary surface. 

(c) Equation (27) evaluated on r = r s . 

5. The system of surface equations described above may be closed with the 
additional information about some of the physical variables evaluated 
on the boundary surface (e.g. the luminosity). 

6. Once the system of surface equations is closed, it may be integrated for 
any particular initial data. 

7. Feeding back the result of integration in the expressions for m and v, 
these two functions are completely determined. 

8. With the input from the point 7 above, and remembering that once 
metric functions are fully specified, field equations become an algebraic 
system of equations for the physical variables; these may be found for 
any piece of matter distribution. 
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4.4 The Surface equations 

As it should be clear from the above, the crucial point in the algorithm is 
the system of surface equations. So, let us specify them now. 

Introducing the dimenssionles variables 

A = r s /m E (0), 

F = 1 - 2M/A, 
M = m s /m s (0), 
Q = u; s , 
(3 = i/m E (0), 

where ms(0) denotes the total initial mass, we obtain the first surface equa- 
tion by evaluating (12) at r = r s . Thus, one gets 

% - ™- < 69 > 

Next, using junction conditions, one obtains from (53), (13) and (16) 
evaluated at r = rs, that 

^ = + n)£, (70) 

with 

E = 47rr|e s , (71) 
24 



where the first and second term on the right of (70) represent the gravitational 
redshift and the Doppler shift corrections, respectively. 

Then, defining the luminosity perceived by an observer at infinity as 

dM 



dp 



we obtain the second surface equation in the form 



f r l(l-F )n + 2LIA. 



(72) 



The third surface equation may be obtained by evaluating at the bound- 
ary surface the conservation law T£ = 0, which reads 

(p + P)(47rr 3 P + m) 



P + 



r(r — 2m) 



4irr(r — 2m) 



m + 



3m 2 



+- (P-P ). 

r-2m 2 r K J 



(73) 



In the case when the effective density is separable, i.e., p = jF(t)H(r); 
equation (73) evaluated at the boundary surface leads to 



dtt 

d(3 



ti 2 



8F_ 

~A 

F 



+ 2F/C(r s ) + 47rp s A(3 - Q 2 ) 



where 



R 



1 — 2m /r ^ r 2 ' 



(74) 



(75) 
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E = E(l + n) (76) 

and 

/C(r E ) = -jt-ln (J- jf* drr 2 H{r)/H{r^j . (77) 

Before analyzing specific models, some interesting conclusions can be ob- 
tained at this level of generality. One of these conclusions concerns the 
condition of bouncing at the surface which, of course, is related to the oc- 
currence of a minimum radius A. According to (69) this requires Q = 0, and 
we have 

w ~ w (78) 

or using (74) 

WO w T 2E 1 

R+^^-r ■ (79) 



dp p s 



inr^A 



Observe that a positive energy flux (E) tends to decrease the radius of the 
sphere, i.e., it favors the compactification of the object, which is easily un- 
derstood. The same happens when R > 0. The opposite effect occurs when 
these quantities have the opposite signs. Now, for a positive energy flux the 
sphere can only bounce at its surface when 

> = 0)>0. 
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According to (79) this requires 

-R(Q = 0) > 0. (80) 

A physical meaning can be associated to this equation as follows. For 
non-radiating, static configuration, R as defined by (75) consists of two parts. 
The first term which represents the hydrodynamical force (see (28)) and the 
second which is of course the gravitational force. The resulting force in the 
sense of increasing r is precisely — R, if this is positive a net outward accel- 
eration occurs and vice-versa. Equation (80) is the natural generalization of 
this result for general non-static configurations. 

5 Models and their numerical implementa- 
tion 

5.1 Effective variables 

Once the profiles of energy density and pressure have been established in the 
satic case via the integration of the corresponding Lane-Emden equations, 
we proceed with the determination of effective variables according to the 
algorithm sketched above. However, as it should be clear such determination 
is not unique. The following possibilities arise: 

1. p — f(t) + h(r) p = g{t) + i(r), where h(r) and i(r) correspond to 
the pressure and total energy density obtained from the integration of 
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the relativistic Lane-Emden equations, in both cases described above. 
However this model has not static limit. 

2. p — f(t)h(r) and p = g(t) + Kpl +1 ^ n , for the case I, where p = 
f(t)ho(r), being ho(r) the baryonic mass density in the static limit; 
p = g(t) + Kp 1+1 / n , for the case II. In both cases K = K(m s ,r s ). 

On what follows we shall consider only the possibility 1 above. 

5.2 Numerical implementation of models 

We have used an standard Runge-Kutta routine to obtain functions h(r), 
ho(r) and i(r) from the integration of relativistic Lane-Emden equations for 
different values of n and a. Integration was performed from £ = until the 
first zero of ip (or ip ). 

Next, for the third surface equation we need to calculate numerically the 
following terms: 



Observe that dk(t)/dt = 0, since h(r s ) = 0. 

For the calculation of (81) and (82) we have adjusted a Chebyshev poly- 
nomial [22] to functions h(r) and i(r). Also, for the calculation of either of 



di{r) 



(81) 



dr 



j r=r s 




(82) 
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these functions in the Chebyshev's nodes or within different interior regions 
we have used the interpolating Lagrange polynomials. 

A standard Runge-Kutta method has also been applied to solve surface 
equations. These three equations are solved as an initial value problem, upon 
specification of A(t = 0), F(t = 0) and one function of u. Specifically we 
choose 



L = e -4(t-5/2) 



2 



where ijir is the mass to be radiated. 

Once the surface equations have been integrated, we proceed to calculate 
the metric functions and their derivatives. For doing so, we need to calculate 
numerically the following expressions: 

[\ 2 h{r)dr, (84) 
Jo 

Ja r(r — 2m) 

and 

Ja dt \ r(r -2m) J K ' 

Where the last expression appears in the equation for the time derivative of 

v given by 

. _ . r d f2(47rr 3 p + m)l f 2(47rr 3 p + m) 1 

ir E (t) <9t { r(r — 2m) J [ r(r — 2m) 



rs(t) 
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For the numerical integration of (85) and (86) it is necessary to calculate 
the integrands on points of the latice defined in the integration of Lane- 
Emden equation, using again the Chebyshev's polynomials and Lagrange 
interpolants. Once the metric functions and their derivatives have been com- 
pletely determined, we use the field equations to obtain algebraically the 
physical variables. 

All along evolution we keep radial dependence obtained from the solution 
of the Lane-Emden equations. This was implemented fitting the h(r) and 
i(r) profiles to the radius's distribution at time t. Thus, the radial coordinate 
is scaled by means of: 



We use as many nodes as interior regions studied. One tipical run takes, 
for one region and one time unit, one a half hour in a 900 MHz. central 
processing unit. 



Although a large number of models has been worked out, we shall present 
here only two for illustration, corresponding to the cases I and II. Both were 
calculated for values: n — 2, a — 0.1, f2(0) = —0.05, with an emission of 




The developed code was paralelized using MPI routines for FORTRAN. 



5.3 Models 
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0.01 of the initial mass. The profiles of physical variables are exhibited in 
figures 3-6. As we increase the emission we arrive at a point where case II 
becomes unphysical before case I. If we increase n, for both cases, the initial 
distribution is less compact. On the contrary, if we increase a the initial 
distribution is more compact. Figure 7 shows the normalized radii evolution 
for both cases, different values of n and a — 0.1. 

6 Conclusions 

We have considered two possible definitions of relativistic polytrope and have 
presented a method to study their evolution. The models represent a gener- 
alization of the static polytrope to the case of evolving and dissipating fluid 
spheres, which in the static limit satisfy a polytropic equation of state. This 
allows for modeling self-gravitating objects, and at the same time brings 
out differences between the two possible definitions of polytropes, considered 
here. We have incorporated dissipation, a fundamental process in the process 
of gravitational collapse, into discussion. It remains, to consider the diffu- 
sion limit, however because of the additional complication associated to the 
necessity of introducing an equation of transport, we have only considered 
here the simplest, streaming out limit. 

Although the examples are presented with the sole purpose of illustrating 



31 



the method (our main goal here being to provide a tool for modeling the 
evolution of relativistic polytropes), some comments on them, are in order. 

The difference between the two definitions of polytrope considered here, 
are clearly exhibited in figures 1-2. To magnify such difeference we present 
the results corresponding to the "ultra-relativistic" case (a — 1). As can 
be seen, for n < 1.365, configurations of case I have smaller radii than those 
corresponding to the case II. This situation reverses for n > 1.365. In general, 
bounded relativistic configurations exist for smaller values of n , than in the 
Newtonian case. 

Figures 3-6 show how differently, both polytropes evolve. As it appears 
from these figures, the case II leads to an stronger collpase. This tendency is 
confirmed by curves c — d of figure 7. Also, curves a — b on this same figure 
show an example of bouncing for n = 2.5. The strongest bouncing of case 
I, further indicates that the equation of state resulting from case I is stiffer 
than the obtained from case II. It is worth mentioning that these differences 
are observed in a large number of models, for a wide range of values of n, a 
and initial data. 
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0.5 1 1.5 2 2.5 

Figure 1: ipo (Case I) and ^ (Case II) as a function of £ for a = 1 and 
different values of n: (a) Case I, n — 0.5; (b) Case II, n = 0.5; (c) Case I, 
n = 1.5; (d) Case II, n = 1.5. 
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Figure 2: ip (Case I) and ip (Case II) as a function of £ for a = la 
different values of n: (a) Case II, n = 2; (b) Case I, n = 2; (c) Case 
n = 2.5; (d) Case I, n = 2.5. 
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Figure 3: Adimensional energy density (Case I to the left; Case II to the 
right) as a function of dimenssionles time for n = 2 and a — 0.1 at different 
regions: (a) r/a = 0.25 (multiplied by 10); (b) r/a = 0.50 (multiplied by 10); 
(c) r/a = 0.75 (multiplied by 10 2 ); (d) r/a = 1.00 (multiplied by 10 4 ). 
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Figure 4: Adimensional pressure (Case I to the left; Case II to the right) as 
a function of dimenssionles time for n = 2 and a — 0.1 at different regions: 
(a) r /a = 0.25 (multiplied by 10 3 ); (b) r/a = 0.50 (multiplied by 10 2 ); (c) 
r/a = 0.75 (multiplied by 10 3 ); (d) r/a = 1.00. 
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Figure 5: Radial velocity (Case I to the left; Case II to the right) as a 
function of dimenssionles time for n = 2 and a — 0.1 at different regions: (a) 
r/a = 0.25; (b) r/a = 0.50; (c) r/a = 0.75; (d) r/a = 1.00. 



40 



. 5 



. 5 




-3 1 1 1 1 1 1 1 1 1 1 1 -3 1 1 1 1 1 1 1 1 1 1 1 

0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 

P P 

Figure 6: Adimensional flux of energy (Case I to the left; Case II to the 
right) as a function of dimenssionles time for n = 2 and a = 0.1 at different 
regions: (a) r/a = 0.25 (multiplied by 10 3 ); (b) r/a = 0.50 (multiplied by 
10 2 ); (c) r/a = 0.75 (multiplied by 10 3 ); (d) r/a = 1.00 (multiplied by 10 4 ). 
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Figure 7: Evolution of the normalized radii for a = 0.1, both cases and 
different values of n: (a) Case I, n — 2.5; (b) Case II, n = 2.5; (c) Case I, 
n = 1.5; (d) Case II, n = 1.5. 
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